1 Etape 1 : Première manipulation du RPG

1.1 Je choisis un point sur la carte : https://www.google.fr/maps

=> clic droit sur la carte + clic gauche sur les coordonnées (copiées dans le presse-papier)

J’indique un rayon en mètres (< 15000) pour sélectionner les parcelles autour du point

Code
# coord_gmaps<- rstudioapi::showPrompt(title = "Collez les coordonnées", message = "coordonnées Gmaps", default = "")
coord_gmaps<-"43.447894436406216, 1.2886163291688764"
lat<-as.numeric(str_split(coord_gmaps,fixed(","),simplify = TRUE)[,1])
lon<-as.numeric(str_split(coord_gmaps,fixed(","),simplify = TRUE)[,2])
# rayon<-as.numeric(rstudioapi::showPrompt(title = "Rayon", message = "Rayon (en m)", default = ""))
rayon<-10000

# création d'une table sf «point» avec les coordonnées saisies
# transformation des coordonnées en syst de proj 2154 (Lambert II - Français) 
point<-data.frame(lon,lat,rayon) %>% 
  st_as_sf(coords = c("lon","lat"),crs = "EPSG:4326") %>%
  mutate(coord_pt_gps=st_as_text(geometry)) %>% 
  st_transform("EPSG:2154") %>% 
  st_sf() %>% clean_names() %>% 
  rename(geom=geometry)

# st_crs(point)
# st_geometry_type(point)
# plot(point)

La table des parcelles agricoles 2021 se trouve sur un serveur PostgreSQL/PostGIS.

1.2 Je me connecte au serveur PostgreSQL avec le mot de passe disponible sur Onyxia/Mes services/PostgreSQL/Read me

Code
# le mot de passe est stocké dans un secret Vault

# postgresql_password <- rstudioapi::askForPassword(prompt = "Entrez le password PostgreSQL")
# postgresql_password <- "1tfawt3nj7fgzo3w7cma"
# Sys.setenv("PASS_POSTGRESQL"="1tfawt3nj7fgzo3w7cma")
  
# Connection à PostgreSQL
 cnx <- dbConnect(PostgreSQL(),
                user = "projet-funathon",
                password = Sys.getenv("PASS_POSTGRESQL"),
               host = "postgresql-758156",
               dbname = "defaultdb",
                port = 5432,
                options="-c search_path=rpg,public") # specify what schema to connect to

1.3 je lance une requête SQL pour sélectionner les parcelles situées autour de mon point dans le rayon choisi.

Code
# suppression de la table «point» si elle existe
dbSendQuery(cnx,"DROP TABLE IF EXISTS rpg.point CASCADE;")
<PostgreSQLResult>
Code
# écriture de la table point dans la base PostGis
write_sf(point, cnx, append = T)

# ajout d'une clé primaire
dbSendQuery(cnx,"ALTER TABLE rpg.point ADD CONSTRAINT point_pkey PRIMARY KEY(coord_pt_gps);")
<PostgreSQLResult>
Code
# ajout d'un index
dbSendQuery(cnx,"CREATE INDEX ON rpg.point USING gist(geom);")
<PostgreSQLResult>
Code
# dbExecute(cnx,"CREATE INDEX ON rpg.point USING gist(geom);")

# 3 - Exécution de la requête de découpage du RPG autour du point sur PostGis  -------

dbSendQuery(cnx,"DROP TABLE IF EXISTS rpg.parc_prox CASCADE;")
<PostgreSQLResult>
Code
dbSendQuery(cnx,"CREATE TABLE rpg.parc_prox AS
    SELECT row_number() OVER () AS row_id, p.coord_pt_gps, p.rayon, r.*  
    FROM rpg.point p, rpg.parcelles r 
    WHERE ST_DWithin(p.geom,r.geom,p.rayon);")
<PostgreSQLResult>
Code
# ajout d'une clé primaire
dbSendQuery(cnx,"ALTER TABLE rpg.parc_prox ADD CONSTRAINT parc_prox_pk PRIMARY KEY(id_parcel);")
<PostgreSQLResult>
Code
# ajout d'un index
dbSendQuery(cnx,"CREATE INDEX parc_prox_geom_idx ON rpg.parc_prox USING gist(geom);")
<PostgreSQLResult>
Code
# 4 - Téléchargement des parcelles proches sous R-------------------------------

# parc_prox<-read_sf(cnx,parc_prox)

parc_prox<-st_read(cnx, query="select * from rpg.parc_prox;")

1.4 J’affiche les parcelles autour de mon point avec une carte interactive leaflet

Code
#plot(st_geometry(parc_prox))
#plot(st_geometry(point),add = T,col = "red")

# Transformation de la projection car leaflet ne connait que le WGS 84
carte_parc_prox<- parc_prox %>% st_transform(4326)

# Marqueur du point
pt_mark<- point %>% st_transform(4326)

# ajout du libellé des cultures
carte_parc_prox_lib<-carte_parc_prox %>% 
  left_join(lib_cult %>% select(-code_groupe_culture),by=c("code_cultu"="code_culture")) 
#création d'un label ad hoc à afficher en surbrillance au passage de la souris sur la carte.
labels<-sprintf("<strong>id_parcel : </strong>%s<br/>
                <strong>Groupe culture : </strong>%s<br/>
                <strong>Culture : </strong>%s<br/>
                <strong>Surface (ha) : </strong>%s<br/>
                <strong>Département : </strong>%s<br/>
                <strong>Commune : </strong>%s<br/>",
                parc_prox$id_parcel,carte_parc_prox_lib$libelle_groupe_culture,
                carte_parc_prox_lib$libelle_culture,parc_prox$surf_parc,
                parc_prox$insee_dep,parc_prox$nom_com) %>% 
  lapply(htmltools::HTML)

# labels

# création d'une palette de couleurs associée au groupe de culture
factpal <- colorFactor("Paired", parc_prox$code_group)

carte_parc_prox_html <- leaflet(carte_parc_prox_lib) %>% 
   # addProviderTiles("Esri.WorldImagery") %>%
   # addTiles() %>% 
  addTiles("http://wxs.ign.fr/essentiels/wmts?REQUEST=GetTile&SERVICE=WMTS&VERSION=1.0.0&STYLE=normal&TILEMATRIXSET=PM&FORMAT=image/jpeg&LAYER=ORTHOIMAGERY.ORTHOPHOTOS&TILEMATRIX={z}&TILEROW={y}&TILECOL={x}") %>%
  addPolygons( #fillColor="white",
               fillColor=~factpal(code_group),
               weight = 2,
               opacity = 1,
               color = "#ffd966",
               dashArray = "3",
               fillOpacity = 0.5,
               highlight = highlightOptions(
                 weight = 5,
                 color = "#A40000",
                 dashArray = "",
                 fillOpacity = 0.0,
                 bringToFront = TRUE),
               label = labels,
               labelOptions = labelOptions(
                 style = list("font-weight" = "normal", padding = "3px 8px"),
                 textsize = "15px",
                 direction = "auto",
                 encoding="UTF-8")) %>% 
  addMarkers(data=pt_mark,~lon,~lat, popup = ~coord_pt_gps, label = ~coord_pt_gps)
Warning in RColorBrewer::brewer.pal(max(3, n), palette): n too large, allowed maximum for palette Paired is 12
Returning the palette you asked for with that many colors
Code
# leaflet(pt_mark) %>% addTiles() %>% addMarkers(~lon,~lat)

carte_parc_prox_html
Code
# saveWidget(widget = carte_parc_prox_html, file = "carte_parc_prox_html.html")

1.5 Quelle est la structure des parcelles agricoles autour de mon point ?

Code
t1 <- parc_prox %>% st_drop_geometry() %>% count(code_group) %>% 
  add_tally(n) %>% 
  mutate(n_pct=round(100*n/nn,1)) %>% 
  select(-nn) %>% rename(n_parcelles=n) %>%
  # adorn_totals() %>% 
cbind(
  # comptage des surfaces
parc_prox %>% st_drop_geometry() %>% count(code_group,wt=surf_parc) %>% 
      add_tally(n) %>% 
      mutate(surf_pct=round(100*n/nn,1)) %>%
      select(-nn) %>%  
      rename(surf_parc_ha=n) %>% select(surf_parc_ha, surf_pct) # %>% 
       # adorn_totals()
) %>% left_join(lib_group_cult,by=c("code_group"="code_groupe_culture")) %>% 
  select(code_group,libelle_groupe_culture,everything()) %>% 
  # arrange(as.numeric(code_group)) %>% 
  arrange(desc(surf_parc_ha)) %>% 
  adorn_totals() %>% 
  mutate(taille_moy_parc=round(surf_parc_ha/n_parcelles,1))
Storing counts in `nn`, as `n` already present in input
Storing counts in `nn`, as `n` already present in input
ℹ Use `name = "new_name"` to pick a new name.
Code
t1  %>% 
  setNames(c("code","groupe de cultures","nombre de parcelles","(%)","surface (ha)","surface (%)","taille moyenne (ha)")) %>% 
  kable(
    format="html",
    caption="<span style='font-size:medium'>Groupes de cultures <strong>locales</strong> par surfaces décroissantes</span>",
    format.args = list(decimal.mark = ",", big.mark = " "),
    booktabs = TRUE) %>%
  kable_styling(font_size = 15) %>% 
  gsub("font-size: initial !important;",
       "font-size: 20pt !important;",.)%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>% 
  row_spec(nrow(t1), bold = T, color = "white", background = "grey")
Structure des cultures au niveau local
code groupe de cultures nombre de parcelles (%) surface (ha) surface (%) taille moyenne (ha)
1 Blé tendre 631 10,9 3 044,52 17,6 4,8
2 Maïs grain et ensilage 368 6,4 2 922,39 16,9 7,9
4 Autres céréales 404 7,0 1 826,24 10,6 4,5
11 Gel (surfaces gelées sans production) 1 437 24,9 1 494,91 8,7 1,0
6 Tournesol 308 5,3 1 430,40 8,3 4,6
18 Prairies permanentes 573 9,9 1 260,24 7,3 2,2
7 Autres oléagineux 184 3,2 1 237,37 7,2 6,7
3 Orge 214 3,7 947,61 5,5 4,4
16 Fourrage 234 4,1 804,31 4,7 3,4
19 Prairies temporaires 302 5,2 802,74 4,6 2,7
8 Protéagineux 116 2,0 538,13 3,1 4,6
5 Colza 106 1,8 468,59 2,7 4,4
28 Divers 754 13,1 184,22 1,1 0,2
15 Légumineuses à grains 34 0,6 147,28 0,9 4,3
25 Légumes ou fleurs 38 0,7 46,95 0,3 1,2
21 Vignes 21 0,4 37,31 0,2 1,8
17 Estives et landes 16 0,3 33,19 0,2 2,1
24 Autres cultures industrielles 23 0,4 32,70 0,2 1,4
20 Vergers 7 0,1 4,76 0,0 0,7
Total - 5 770 100,0 17 263,86 100,1 3,0
Code
# rm(t1)

1.6 Comparaison avec la répartition des cultures au niveau départemental et national

Code
# jointure spatiale avec la couche département pour récupérer le département où tombe le point

# ouvrir la couche des communes (à convertir en Lambert II epsg 2154)
dep <- s3read_using(
    FUN = sf::read_sf,
    layer = "departement",
    object = "2023/sujet2/diffusion/ign/adminexpress_cog_simpl_000_2023.gpkg",
    bucket = "projet-funathon",
    opts = list("region" = "")) %>% 
  st_transform(2154)

# jointure
df<-point %>% st_join(dep) %>% st_drop_geometry() %>% select(insee_dep)
dep_pt<-df[1,1]
rm(df)

# sinon sélection du département regroupant la plus grande surface agricole
# cas où le cercle recouvre plusieurs départements
# df<-parc_prox %>% st_drop_geometry() %>% 
#   count(insee_dep,wt=surf_parc) %>% 
#   arrange(desc(n))
# dep_pt<-df[1,1]
# rm(df)

# calcul des % surfaces autour du point
stat_pt <- parc_prox %>% st_drop_geometry() %>% 
  count(code_group,wt=surf_parc) %>% add_tally(n) %>% 
  mutate(pct_surf_local=round(100*n/nn,1)) %>%
  select(code_group, pct_surf_local) 
Storing counts in `nn`, as `n` already present in input
ℹ Use `name = "new_name"` to pick a new name.
Code
# récup des % surfaces départementales
stat_dep_pt<-s3read_using(
  FUN=readr::read_rds,
  object = "2023/sujet2/diffusion/resultats/stat_group_cult_by_dep.rds",
  bucket = "projet-funathon",
  opts = list("region" = "")) %>% 
  filter(insee_dep %in% dep_pt) %>% 
  select(insee_dep,code_group,libelle_groupe_culture,pct_surf) %>% 
  rename(pct_surf_dep = pct_surf)

# récup des % surfaces nationales
stat_fm<-s3read_using(
  FUN=readr::read_csv,
  object = "2023/sujet2/diffusion/resultats/stat_group_cult_fm.csv",
  col_types = cols(code_group = col_character()),
  bucket = "projet-funathon",
  opts = list("region" = "")) %>% 
 select(code_group,libelle_groupe_culture,pct_surf) %>% 
  rename(pct_surf_fm = pct_surf)

# appariement des stas locales, départementales, nationales
stat_compar<-stat_fm %>% 
  left_join(stat_dep_pt %>% select(code_group, pct_surf_dep),by="code_group") %>% 
  left_join(stat_pt,by="code_group") %>% 
  select(libelle_groupe_culture,pct_surf_local,pct_surf_dep,pct_surf_fm) %>% 
  arrange(desc(pct_surf_local)) %>% adorn_totals() 

stat_compar %>% 
  setNames(c("Groupe de cultures","surf. locales (%)","surf. départ. (%)","surf. France m. (%)")) %>% kable(
    format="html",
    caption="<span style='font-size:medium'>Comparaison des surfaces locales, départemenales et nationales</span>",
    format.args = list(decimal.mark = ",", big.mark = " "),
    booktabs = TRUE) %>%
  kable_styling(font_size = 15) %>% 
  gsub("font-size: initial !important;",
       "font-size: 20pt !important;",.)%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>% 
  row_spec(nrow(stat_compar), bold = T, color = "white", background = "grey")
Structure des cultures
Groupe de cultures surf. locales (%) surf. départ. (%) surf. France m. (%)
Blé tendre 17,6 14,3 17,7
Maïs grain et ensilage 16,9 8,6 10,0
Autres céréales 10,6 12,4 3,9
Gel (surfaces gelées sans production) 8,7 4,4 1,6
Tournesol 8,3 12,6 2,5
Prairies permanentes 7,3 17,9 27,7
Autres oléagineux 7,2 3,1 0,7
Orge 5,5 3,5 6,2
Fourrage 4,7 5,9 3,7
Prairies temporaires 4,6 4,2 5,3
Protéagineux 3,1 1,8 1,2
Colza 2,7 1,8 3,5
Divers 1,1 1,3 1,3
Légumineuses à grains 0,9 0,3 0,2
Légumes ou fleurs 0,3 0,3 1,5
Estives et landes 0,2 6,9 8,1
Vignes 0,2 0,4 2,2
Autres cultures industrielles 0,2 0,1 1,8
Vergers 0,0 0,1 0,4
Plantes à fibres NA 0,0 0,4
Riz NA NA 0,0
Fruits à coque NA 0,0 0,2
Oliviers NA 0,0 0,0
Total 100,1 99,9 100,1

1.7 Graphique de comparaison des cultures au niveau local, départemental et national

Code
# je sélectionne les 10 groupes de cultures les plus répandus au niveau local 
tab<-stat_compar %>% filter(libelle_groupe_culture!="Total") %>% slice_head(n=10) %>% 
  rename(local=pct_surf_local, departement=pct_surf_dep, france=pct_surf_fm)

# je transpose la table pour rassembler toutes les valeurs dans une seule variable value (ggplot oblige)
tab_piv<-tab %>% pivot_longer(!libelle_groupe_culture) %>% rename(secteur=name) 

# je mets à 0 les valeurs manquantes
#tab_piv<-tab_piv %>% coalesce(value,0L)  
tab_piv[is.na(tab_piv)] <- 0


# levels(as.factor(tab_piv$secteur))
# je réordonne les secteurs dans le "bon" ordre, avec factor  
tab_piv$secteur<-factor(tab_piv$secteur,levels=c("france","departement","local"))
tab_piv<-tab_piv %>% arrange(desc(secteur), desc(value))

# je réordonne les cultures par surface décroissante au niveau local, avec factor
x <- tab_piv %>% filter(secteur == "local") %>% arrange(value) %>% select(libelle_groupe_culture)
y <- pull(x, libelle_groupe_culture)

tab_piv$libelle_groupe_culture <- factor(tab_piv$libelle_groupe_culture, levels = y)

# ggplot => problème de tri : affichage des prairies avant le maïs ?    
p<-ggplot(tab_piv, aes(x =libelle_groupe_culture,
                       y = value, 
                       fill=factor(secteur, levels=c("france","departement","local")))) + 
  geom_col(position = "dodge") +
  labs(title="Surfaces comparées des 10 principales cultures locales, en %", x="Culture", y = "%", fill = "Secteur") +
  theme_classic()

# je flippe le graphique pour avoir des barres horizontales  
p+coord_flip()

différences de structure locale, departementale, nationale

1.8 Je teste un graphique par secteur avec facet_wrap

Code
# je sélectionne les 10 groupes de cultures les plus répandus au niveau local 
tab<-stat_compar %>% filter(libelle_groupe_culture!="Total") %>% slice_head(n=10) %>% 
  rename(local=pct_surf_local, departement=pct_surf_dep, france=pct_surf_fm)

# je transpose la table pour rassembler toutes les valeurs dans une seule variable value (ggplot oblige)
tab_piv<-tab %>% pivot_longer(!libelle_groupe_culture) %>% rename(secteur=name) 

# je mets à 0 les valeurs manquantes
#tab_piv<-tab_piv %>% coalesce(value,0L)  
tab_piv[is.na(tab_piv)] <- 0

# ggplot => problème de tri : affichage des prairies avant le maïs ?    
ggplot(tab_piv, aes(x = reorder(libelle_groupe_culture,+value), 
                       y = value)) + 
  geom_col(fill = "lightblue", colour = "black",position = "dodge") +
  labs(title="Surface par culture", x="Culture", y = "%", fill = "Secteur") +
  geom_text(aes(label = value), hjust = -0.3, size = 8/.pt, colour = "black") +
  theme_classic() + coord_flip() + 
  facet_wrap(~secteur,nrow=3,scales='free')

différences de structure locale, departementale, nationale

1.9 Pour aller plus loin, comparer la taille moyenne des parcelles selon les secteurs